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A  MODIFIED  SUBSONIC  AERODYNAMIC  SUBROUTINE  FOR 
THE  VIBRA  NUCLEAR  VULNERABILITY  CODE 

1.  INTRODUCTION. 

The  aerodynamic  subroutine  in  VIBRA-4  (Reference  1)  for  lifting 
surfaces  of  aircraft  flying  at  subsonic  speeds  was  devised  from  empirical 
results.  For  all  practical  purposes,  this  subroutine  is  essentially  the 
same  as  that  in  the  earlier  VIBRA  versions.  Prepared  nearly  twenty 
years  ago,  the  original  version  (Reference  2)  was  tailored  specifically 
for  configurations  of  interest  at  that  time.  The  data  utilized  in 
Reference  2  consist  of  "average"  values  which  are  deemed  acceptable 
whenever  the  Mach  number  (M)  is  around  0.8,  the  aspect  ratio  (AR)  is 
between  4  and  6,  and  the  quarter-chord  sweep  (A)  is  say  between  35  and 
45  degrees.  Some  current  aeronautical  systems  have  higher  aspect  ratios, 
lower  sweeps,  and  the  flight  speeds  are  such  that  M  falls  between  0.5 
and  0.75. 

In  what  follows,  certain  practical  modifications  are  effected 
which  allow  the  subsonic  aerodynamic  subroutine  to  be  applicable  over 
broader  ranges  of  Mach  number,  aspect  ratio,  and  sweep;  and,  more 
specifically,  over  the  ranges  0.5  _<  M  <_  0.8,  3  _<  AR  £  12,  and  0  <_  A  _<  60 
degrees.  The  latest  VIBRA  version,  designated  as  VIBRA-8,  differs  from 
VIBRA-4  in  that  it  includes  the  modified  and  improved  subsonic  aero¬ 
dynamic  subroutine.  Section  II  of  the  original  VIBRA  report  (Reference  2), 
which  documents  the  subsonic  formulation  for  all  versions  through 
VIBRA-4,  is  followed  in  the  discussions  below.  The  intent  of  this  short 
report  Is  to  justify  and  to  document  the  effected  modifications  and 
additions. 

2.  AERODYNAMIC  CONSIDERATIONS. 

There  are  three  areas  which  need  to  be  examined: 

1.  Potential  flow  results  for  the  ratio  of  two-dimensional 
normal  force  coefficient  to  angle  of  attack  as  a  function 
of  aerodynamic  time. 
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2-1 


2.  Two-dimensional,  non-potential  flow,  large  angle  of 
attack  effects. 

3.  Three-dimensional  effects. 


Ratio  of  Two-Dimensional  Normal  Force  Coefficient  to  Angle 
of  Attack  as  a  Function  of  Aerodynamic  Time. 

Consider  first  the  function  f(s). 


f(s) 


cn(a,s) 

a 


(1) 


which  is  the  ratio  of  the  normal  force  coefficient  to  angle  of  attack 
according  to  linearized  potential  flow  theory  for  a  two-dimensional 
unswept  wing  subjected  to  a  step  change  in  angle  of  attack,  i.e.,  the 
indicial  problem.  This  functionality  with  respect  to  the  dimensionless 
time  parameter  s  is  given  by  Eq.  (128)  and  Figure  9  in  Reference  2,  with 
the  units  of  f (s)  being  per  degree.  The  value  of  f (s)  according  to 
Eq.  (128)  should  be  related  to  the  so-called  Wagner  function,  <t>c,  as  given 
on  Page  350  of  Reference  3.  In  fact,  it  can  be  easily  shown  that 


f(s)  *  2it$c(s)  (2a) 

if  f (s)  is  expressed  in  units  of  per  radian,  and 

2 

f(s)  =  fo*c(s)  (2b) 

if  f(s)  is  expressed  In  units  of  per  degree. 

Invariably,  the  function  <j>c  is  approximated  in  practical 
applications  through  the  expression 

...  ,  .  -26,s  ,  -26-8  .  .  -28-s 

<t>c(s)  *  bQ  +  b^e  1  +  b2e  2  +  b3e  3  (3) 

where  b^,...b3>  are  parameters  depending  on  Mach  number  only. 

This  expression  differs  slightly  from  the  corresponding  one  in  Reference  3. 
The  adjustment  here  is  due  to  the  fact  that  in  Reference  3  the  parameter 
s*Ut/c,  whereas  in  Reference  2  and  here  s*2Ut/c,  where  U  is  the  wind 
speed,  t  is  time,  and  c  is  the  chord. 


From  steady-state  considerations  for  an  unswept  two-dimensional 
(2D)  airfoil  in  subsonic  compressible  flow. 


c  (s=<=) 

n 

a 


f  (S=») 


(4) 


It  is  also  known  that  0^,  $3  are  always  positive;  and  thus  terms 

involving  them  -*•  0  as  s  -*  «>.  From  Eq.  (4)  and  Eq.  (2a), 


4>c(s-+“) 


(5) 


From  acoustic  theory  (or  linearized  piston  theory) ,  the  condition  that 
the  loading  at  t=s=0  depends  on  the  instantaneous  angle  of  attack,  and 
more  precisely. 


cn(s=0)  4 

— - -  =  f  (s=0)  =  —  (per  radian) 

leads  to  the  result  that 

(b0  +  bl  +  b2  +  b3jmust  eclual  to 

or  equivalently. 


b 


3 


(6) 


except  at  M=0  .  The  entries  in' Table  6-1  of  Reference  3  are  In  accord 
with  the  above  results. 


To  avoid  confusion  in  later  derivations,  quantities  obtained 
from  potential  flow  theory  are  so  designated  by  a  superscript  P.  Thus, 
when  expressed  in  units  of  per  degree. 


fP(s) 


(s,a-*0)j 
90  (b0  +  bl 


“20. s  .  -20  s  -20,s 

e  1  +  b^e  l  +  b^e  2  } 


(7) 


When  M*0,  one  has  to  deal  with  an  impulse  due  to  apparent  mass  effects. 
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It  should  be  noted  that  f(s)  as  given  by  Equation  (128)  of  the  VIBRA 

p 

report  is  really  f  (s);  and,  as  far  as  can  be  determined,  it  was  obtained 
from  data  for  a  single  Mach  number,  namely  M=0.8.  The  first  task  at 
hand  is  to  give  formulas  for  b^,...,  6^  as  functions  of  M,  so  that  the 
more  general  equation  (i.e.,  Eq.  (7)  above  instead  of  Eq.  (128))  may 

p 

be  applied  to  represent  f  (s)  for  a  range  of  Mach  numbers,  and  not  be 
limited  to  the  single  M=0.8.  Using  the  numerical  values  from  Table  6-1 
of  Reference  3  and  corresponding  values  for  M=0.8  from  Reference  4, 
curves  were  fitted  for  each  of  the  parameters  b^,...,8^.  The  curve  fits 
resulted  in  the  following  expressions  for  the  range  0.2  £  M  £  0.8 

1 

/l-M2 


b 

o 


b 


1 


b 


2 


b 

8 


3 

1 


8 


2 


8 


3 


-0.165  -  0.482M 

1.044  -  6.933333M  +  11. 4M2  -  6.666667M3 

-0.249  -  0.344  (0.5  -M)2 

36.806  -  173.393M  +  266. 4M2  -  135.668M3 


0.0754  -  0.1196  (0.5  -  M)2 

0.4029  -  1.50583M  +  2.36M2  -  1.31667M3 

0.300  +  0.288M2 

-13.293  +  62.94667M  -  94.15M2  +  45.83333M3 
0.945M-1 

58.96  -  257. 59M  +  380. 4M2  -  187. OM3 


0.2  £  M  <  0.5 
0.5  £  M  £  0.8 
0.2  £  M  <  0.5 
0.5  £  M  £  0.8 

0.2  £  M  <  0.5 
0.5  £  M  £  0.8 
0.2  £  M  £  0.5 
0.5  £  M  £  0.8 
0.2  £  M  £  0.5 
0.5  <  M  <  0.8 


Through  extensive  calculations,  it  has  been  determined  that  when  Eq.  (7) 
is  applied  using  the  curve  fitting  expressions  (Eq.  (8))  considerable 

p 

improvements  are  achieved  in  estimating  f  (s) ;  this  is  true  even  at 
the  single  Mach  number  (*0.8)  for  which  Eq.  (128)  of  Reference  2  was 
offered  earlier. 


(8) 
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/ 


made: 


2-2 


Based  on  the  above  brief  discussions,  the  following  change  was 


Change  1: 

(a)  Omit  Equation  (128)  of  Reference  2. 

(b)  To  replace  it,  compute  first  b^,...,^  according  to 
Eq.  (8)  for  the  actual  M. 

(c)  Use  the  b^,...,^  coefficients  thus  computed  to  find 

P.  .  it 2  (  -28,s  -28  s  ,  ,  -28-s( 

f  (S)  =  90  b0  +  V  1  +  b2e  2  +  V  3  | 


(9) 


which  is  in  units  of  per  degree  and  ready  to  replace 
Eq.  (128).  This  eliminates  the  use  of  the  same  and  less 
accurate  formula  for  all  subsonic  M's. 

Two  Dimensional,  Non-Potential  Flow,  Large  Angle  of  Attack 
Effects . 


Refer  once  more  to  Figure  9  of  Reference  2.  If  the  loadings 

p 

were  linear  with  angle  of  attack,  the  f  (s)  as  modified  above  would  suffice 

for  the  two-dimensional  case.  However,  steady-state  measurements  on 
cn 

airfoils  show  than  —  (s-*>°,  a/0)  depends  on  a.  This  is  evident  from 
the  asymptotes  (i.e.,  as  s  -+  large  values)  of  the  non-potential  curves 
in  Figure  9.  The  non-potential  conditions  will  henceforth  be  referred 
to  as  the  NP-conditions;  and  to  stress  them,  when  necessary  to  do  so  for 
the  sake  of  clarity,  related  quantities  will  be  superscripted  with  NP. 


As  seen  in  the  figure,  the  potential  flow  conditions  prevail 

up  to  times  s=s,  (a);  thereafter,  the  curves  decay  towards  their  steady- 

°  NP 

state  values  as  s-*=°.  Equation  129  of  Reference  2  gives  f  (s)  in  the  form 


fNP(s) 


A  +  BAe_CAS 
A  A 


(10) 


where  A.,  B.,  and  C„  are  functions  of  u.  Table  I  of  the  same  reference 
A  A  A 

gives  the  variations  of  A  ,  B  ,  and  C  with  u.  Presumably,  these  have 

AAA 

been  obtained  from  shock  tube  data  and  correspond  to  M=0.8.  This  f  (s) 
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is  used  for  s  >  s.  in  previous  VIBRA  versions.  Equation  (10)  needs  to 
t> 

be  reexamined.  Let  it  be  cast  into  the  form 


cNP,  \ 

f  (s)  *  A, 


K") 


A  (a-0) 
A 


'  AA(a#0) 

A  (a=0) 
i  A 


The  first  factor  on  the  right  hand  side,  i.e.,  A^(a=0),  is  of  special 

irAA(o=0) 

significance.  The  quantity  — jgg -  which  is  A^(a-O)  of  Table  I, 

p 

but  in  units  of  per  degree  rather  than  per  radian,  is  really  (c^/a)  , 
steady-state.  Based  on  potential  flow  derivations  in  Reference  3, 

/  Cn  \  _  *  _  *2  1 


ft) 


_  TT  -  ,  TT _ 

steady— state  180  1TD0  90 


/l-M2 


which  is  equal  to  0.1828  for  M=0.8.  The  entry  AA(a=0)  in  Table  I  gives 
for  the  same  quantity  the  close  value  of  0.1837.  Thus,  an  analytical 
expression  has  been  established  for  AA(a=0,M) , 


AA(a=0,  M)  = 


(per  degree) 


which  is  no  longer  restricted  to  Mach  numbers  near  0.8. 


Attention  is  next  focused  on  the  second  term  which  is  the  ratio 


r(a)  =  r 


Aa(“#0) 

AA(a*0) 


(13) 


NP 

Note  that  f  -  AA(af*0)  *  (AA(a*0))  r.  From  Table  I,  Reference  2, 

the  following  values  are  obtained  for  r(a,  M=0.8): 


«(deg) 


0 

5 

10 

15 

30 

45 

60 


1.0 

0.6984 

0.4573 

0.3059 

0.2559 

0.2243 

0.1998 
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A  test  was  made  to  find  out  if  this  "distribution  function"  r  fits 
experimental  data  at  other  Mach  numbers.  Had  this  test  shown  a  reason¬ 


able  fit,  one  could  then  use  the  above  r-entries  along  with  A  (a=0,  M) 

NP  A 

as  given  by  Eq.  (12)  to  obtain  f  (s-*») .  Data  for  the  symmetric 

64A006  and  64A010  airfoils  (6  and  10  percent  thick,  respectively) 

presented  in  Reference  5  were  found  suitable  for  this  purpose.  If  one 

plots  the  experimental  data  for  r  =  /cWa)ct^Q  versus  a>  an<j  corresponding 

(Cn/a)a=0 


to  M=0.3,  0.5,  0.7,  and  0.8,  the  various  M-curves  do  not  fall  close  to 
each  other.  However,  if  one  plots  this  ratio  r  as  a  function  of  a ,  where 
is  the  "equivalent  angle  of  attack  at  M=0.8" 

0 . 6a 

a  =  -  , 

e 


then  all  the  data  points  fall  close  to  an  "average  curve".  Since  at  M=0.8, 
ae=a,  this  average  curve  should  correspond  to  r(a=«e,  M=0.8)  and  more  or 
less  duplicate  the  curve  one  has  from  the  r(a,  M=0. 8) -tabulation  above, 
i.e.,  from  Table  I  of  Reference  2.  There  are  some  differences  however. 

The  average  curve  resulting  from  the  steady-state  data  of  Reference  5 
is  better  defined  and  is  therefore  used  here. 


The  above  observations  indicate  that  the  dependency  of  r  on 
two  parameters  a  and  M  can  be  reduced  to  a  dependency  on  the  single 
parameter  a^.  Thus, 


r(a,M)  =  r(ae) 
from  which  it  follows  that 


AA(a,M) 


AA(a=0,  M)  r  (a,M) 


=  A  (u=0,  M)  r  (a  ) 

A  A  £ 

,  ,  /  1 

r(“^  75? 


(14) 


-3 


The  average  curve  representing  r(afi)  vs  and  referred  to  earlier  is 
described  by  the  following  table: 

ag(deg)  rfa^)  from  "Average  Curve" 


0 

1.0 

±5 

0.875 

±10 

0.440 

±15 

0.310 

±30 

0.256 

±45 

0.224 

±90 

0.151 

Linear  interpolation  may  be  used  for  cte's  between  the  tabulated  values. 
This  table  should  be  compared  with  the  preceding  r(a,  M=0.8)-tabulation. 

Consider  next  the  third  factor  on  the  right  hand  side  of 

Eq.  (11).  One  way  of  handling  this  "decay"  is  to  assume  that  the 

NP 

non-potential  value  f  breaks  from  the  potential  value  at  the  breakpoint 

s=s^(a) ,  where  s^'s  are  as  in  the  previous  VIBRA  versions,  and  that  the 

decay  is  such  that  following  30  chord  lengths  of  travel,  i.e.,  at 

NP 

s=s,  +  30,  the  decay  to  the  steady-state  value  f  (a,  s-*»)  is  99% 
b 

complete.  However,  this  is  at  most  a  rough  approximation.  One  would 

NP 

do  just  as  well  to  retain  the  previous  f  -format,  i.e., 


f 


NP 


A  +  Be  CAS 
A  A 


using  B  =  B  (a,  M=0.8)  and  C.  =  C.(ot,  M=0.8)  as  given  in  Table  I  of 

A  A  A  A 

Reference  2.  (Of  course,  it  may  be  argued  that  perhaps  a  better  approach 
would  be  to  say  B^  =  BA(ag)  and  CA  =  CA(ae>,  in  analogy  with  the  change 
associated  with  AA(a,  M) ,  i.e.,  Eq.  (14).  But  there  is  no  applicable  data 
to  support  this  hypothesis.) 


Based  on  the  above  discussions,  the  following  change  was  made: 


Change  2: 

(a)  Use  Eq.  (129)  of  Reference  2,  but  with  AA  (per  radian) 

from  Table  I  of  the  same  reference  replaced  by 
2tt 

A.(a,M)  =  r(a  )  -  (units  here  are  per  radian) 

A  6  /Ihm? 

to  obtain  fN^. 
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(b)  Proceed  In  exactly  the  same  manner  as  before,  but  with  the 
P  NP 

modified  f  (Change  1)  and  f  (Change  2a),  to  obtain 
c 

n  ,  . 

—  (s,a). 

a 

2-3  Three-Dimensional  Effects. 

Attention  is  next  turned  to  methods  for  accounting  for  three- 
dimensional  (3D)  and  sweep  effects.  Consider  first  the  steady-state 
case  which  corresponds  to  s-**. 

CN 

Figure  10A  of  Reference  2  gives  the  ratio  —  (a,s-x»)  as  a  function 

c 

n 

of  a  for  a  "representative  case",  namely  that  for  M  =  0.8,  A  =  35  to  45  deg, 

AR  =  4  to  6.  Of  course,  M,  AR,  and  A  are  expected  to  affect  this  curve. 

Cn 

One  can  see  this  by  examining  the  starting  point,  i.e.,  —  (ot=0,  s->°°) . 

n 

Experimental  data  from  a  NASA  Memorandum  (Reference  6)  have  been  used  to 
confirm  the  well-known  expression  given  as  Eq.  6-35  in  Reference  3  to  be 
quite  accurate.  This  expression  (for  the  normal  force  slope)  is: 


a  cosA 


AR  /l-M2 


AR  /l-M2 


vi+e 


a  cosA 
o  e 


ttAR  /l-M2  i 


where 


A  =  tan 
e 


-1  /  tan  A  \ 
\  /l-M2  / 


(15a) 


The  2D  incompressible  lift  curve  slope  a  may  be  taken  equal  to  2tt,  and  its 
compressible  counterpart  c  =  o  =  2tt  . 


Since 


—  (a=0 ,  s-*»)  = 

n 


one  has  from  Eq.  (15) 


/l-M2'  /l-M2 


—  (a**0,  8-*»)  *  cosA 

c  e 

n 


AR  /l-M2 


AR  /l-M2 


2cosA 


r : 


2cosA 


./ 


N 

Calculations  for  —  (a*0,  b-*»)  according  to  Eq.  (16)  show  large  variations 
n 

with  respect  to  K,  A,  and  AR.  The  single  representative  value  from 
Figure  10A  of  Reference  2  is  around  0.445  (corresponding  to  M*0 . 8 , 

A  *  35  deg,  AR  -  5.5).  But  this  ratio  is  substantially  higher  in  cruise 
missile  applications  where  one  usually  has  lower  M's,  lower  A's,  and 
higher  AR's.  This  indicates  that  Figure  10A  must  be  adjusted  to  reflect 
its  dependency  on  M,  A,  and  AR. 

CN 

Although  the  steady-state  values  of  —  at  a»0  for  different 

n 

(M,  A,  AR)  combinations  can  be  estimated  adequately  by  the  use  of  Eq.  (16), 
the  question  still  remains  as  to  how  to  predict  the  variation  of 

s 

—  (a;  M,  A,  AR)  versus  a.  After  much  exploration,  the  following  approxi- 
n 

mation  was  found  to  be  in  reasonable  agreement  with  data  used  in  the 
earlier  V1BRA  work  and  with  related  data  from  the  NASA  Memorandum 
(Reference  6) : 


(a,  s-K»;  M,A,  AR)  =  £  [1.7-0. 7  cos  10a  ] 

a  e 

*  5  [2.05-0.35  cos  10a  ] 

e 

-  5  [1.7] 

where  a  ■  0 . 6a/ /l-Mz  as  before,  and 

6  CN 

5  »  -2.  (a=0,  s-*»;  M,  A,  AR) 

n 

which  is  exactly  the  quantity  obtained  from  Eq.  (16). 


0  <  ae  <  18  deg 
18<_  ag  £  36  deg 
ae  1  36  deg 


Up  to  this  point,  one  can  get  the  2D  cn(a,  s-*» )  and  the  3D 

s-*»);  but  the  C^-value  is  for  the  total  lifting  surface.  The  next 
question  to  be  addressed  is  how  to  distribute  the  load  spanwise  such  that 
the  total  loading  will  correspond  to  C^.  A  spanwise  distribution,  denoted 
by  g(n)  (n  being  a  dimensionless  spanwise  coordinate),  is  needed,  such 
that  the  local  loading  per  unit  span  is  /g(n)CN 


V“’ 


(17) 
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I 


In  incompressible  flow,  the  sectional  lift  Jt(n)  is  elliptic  if 
the  chord  varies  elliptlcally  along  the  span.  If  one  sets 


*<n) 


qc(n)  c^  a 
a 


and  c(n) 
at  M*0,  one  would  want  c 


cq  /l-nz ,  where  n  =  |2y/§|,  8  =  span,  cq 


root  chord,  then 


=  constant  =  k  ,  i.e.,  independent  of  n. 
a 


Thus  for  an  elliptic  plan  wing  at  M=0. 

Hn) 


qcQ  (ko/l-n^)a 


For  non-elliptic  wings  at  M=0,  let 


*(n)  =  qc 


k  /l-n^ 
o 

f(n) 


c(n) 


(18) 


where  f(n)  is  to  be  chosen  so  that  if  the  c(n)  becomes  elliptic  f(n) 
cancels  the  (c(n)/co)-variation  to  give  the  A-n2  -type  spanwise  dis¬ 
tribution.  An  elliptic  wing  may  be  thought  of  as  approximately  a  tapered 
wing  with  taper  ratio  (TR)  equal  to 


(TR) 


—  -i 
2 


=  0.571 


because  its  area  is 


tic  s 
o 


c  and  tip  chord  ck  is 
o  t 

I  c+c 


m- 


Sc 


and  the  area  of  a  tapered  wing  with  root  chord 


1  +  ((TR)  -l)n  =  l-0.429n 
e 


(19) 


A  suitable  expression  for  f(n)  is  then 
f(n) 

resulting  in 

&(n)  =  ac  1  k  /l-n^  ^(TR-ljq  I 

qco^ko/ln  1_o.429n  f 

Equation  (19)  does  not  show  the  dependence  of  lift  distribution 

on  the  effective  sweep  parameter  Ag  =  tan  *  (tan  A/ A-m2  )  which  is  known 

to  exist  (e.g..  Reference  7).  (The  Ag-parameter  is  the  same  as  the 

commonly  used  A  -parameter  found  in  the  same  reference.)  With  rra  sweep, 

p 

the  sought  distribution  g(n)  is  insensitive  to  M;  but  the  magnitude  of 
Jt(n)  is  dependent  on  kQ  which  in  turn  depends  on  M  through  the  factor 
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/ 


/ 


1//1-M2.  Let  Afi  be  given  in  degrees.  Instead  of  using  Eq.  (19),  a 
modified  version  would  be  more  appropriate.  Assuming  the  A^-variation 
is  a  linear  one,  let 


t(n)  -  qc0  jk/IV  (i  +  »n2  ^f)J« 


(20) 


where  "a"  is  a  suitable  constant  to  be  determined  through  trial  examples 
Let  further 

.1  / 

,o  ^  ^l-67429n  li  +  ar2 

Then  the  total  lift  L  is 

3/2 


l1  **  [m?  &)) 


dn 


2/o 


2  r  Ky)dy 
o 


s  (q)c  k  la 
^  o  o 


(21) 


c  +c 
o  t 


2c 

— 77—  k  aql 
Co+Ct  ° 

2k 


9  S, 


wing  ~  1+TR 


C,  q  S  .  a 
L  wing 
a 


from  which  one  obtains 


2k 


N  L  1+TR 

a  a 


An  important  parameter,  which  is  also  a  spanwise  distribution 
function  and  will  be  used  shortly,  is 


p(n)  = 


V 

CLCav 


C£  C 
a 

CT  c 
L  av 
a 


where  c  is  the  average  chord.  From  the  preceding  equations,  it  turns 
av 

out  that 


P(n)  = 


c*  c 
a 

Cav 

a 


^  1WIR-1W  1  (22) 
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1  2 

By  expanding  j_o  429n  “  1  +  0.429n  +  (0.429n)  ...»  and  performing 

the  Indicated  integration  in  Eq.  (21),  one  has  also 

I  *  (  0.785398  +  0.196350a  -^)+  (TR-0.571)  ^0.453803  +  0.197011a 

(23) 

The  simple  approximation  Eq.  (22)  is  the  basis  from  which  the 
distribution  function  is  derived.  To  determine  whether  it  is  a  reasonable 
one,  its  results  were  compared  with  those  from  the  extensive  parametric 
study  reported  in  Reference  7.  (This  parametric  study  is  based  on  steady- 
state  linearized  theory  and  covers  wide  ranges  of  M(<1),  A,  AR,  and  taper 
ratio  TR.)  Choosing  the  constant  a  to  be  equal  to  1,  it  was  found  that 
Eq.  (22)  is  indeed  an  adequate  representation  for  (M,  A,  AR,  TR) -combinations 
appropriate  for  a  variety  of  current  subsonic  aircraft.  It  was  therefore 
adopted. 

With  some  algebraic  manipulations,  it  can  be  shown  that  g(n)  is 
related  to  p(n) .  From  that  relation,  it  turns  out  that 


g(n) 


1+TR  /l-n2 
2  l-0.429n 


(l  +  n2 


1 

I 


(24) 


where  I  is  according  to  Eq.  (23)  with  a*l 


For  sectional  lift  coefficients  on  a  3D-ving, 
one  can  use 


e 

n 


fin(n)’ 


e 

^  (a»0,  8-+“) 
n 


~  (o*0,  s-*»)  g(n) 
n 


If  one  assumes  that  g(n)  holds  at  all  a's  (since  nothing  better  and  simple 
enough  is  available), 

Cn  CN 

~  (a,  s-«>)  -  —  (a,  s-+«)  g(n)  (25) 

n  n 


The  distribution  of  the  aerodynamic  loading  along  the  span  as 
indicated  by  Eq.  (25)  is  for  8-*°°,  i.e.,  the  steady  state  case.  At  the 


ocher  extreme  of  the  time  scale,  i.e.,  t«s“0,  three-dimensional  effects 
should  be  negligible,  i.e.,  the  flow  is  really  two-dimensional,  and  thus 


e 

^  (a,  s -0)  -  1  (26) 

n 

In  Reference  2,  the  following  procedure  is  suggested  for 
intermediate  times  s: 


C  (a’S) 


1  -  i1  -  -r  (°.  s-~> )  ^ 

'  n  ' 

CN  .  . 

—  (a,  s-**>) 


0  <  s  <  1 


1  <  s  £  16 

s  >  16 


(27) 


(See  Figure  10B  of  Reference  2.)  With  the  introduction  of  the  spanwise 
distribution,  the  following  alternative  to  Eq.  (27)  is  recommended  to 

g 

describe  the  variation  of  _n  with  s: 

c 

n 


e 

(a,  s)  -  1  0  <_  s  <_  1 

n 

-  !  -(l  -«-»•“»-»)(!-  ^  (,.  s~>)  s>! 

e 

Of  course,  —  (a,  s)  must  be  multiplied  by  the  spanwise  strip  width  to 
n 

obtain  the  suitable  coefficient  for  the  strip  i.  The  rest  of  the  VIBRA 
procedure  remains  the  same. 

Based  on  the  above  discussions,  the  following  final  change  and 
addition  was  made: 

Change  3 

(a)  Find  according  to  Eq.  (15a) 

CN 

■  —  (a*0,  s-**>)  according  to  Eq.  (16). 
n 


(b)  Find  £ 


VM 

(c)  Find  —  (a,  s-*«)  according  to  Eq.  (17). 

n 

(d)  Find  I  and  g(n)  for  desired  n's  from  Eq.  (23)  with  a  =  1 
and  Eq.  (24).  Here  n  *  1 2y / g j . 

e 

(e)  Find  —  (a,  s-*=°)  according  to  Eq.  (25). 

c 

n 

This  quantity  will  be  a  function  of  n»  unlike  its  counter¬ 
part  in  earlier  VIBRA  versions. 

C 

(f)  Find,  using  the  results  from  (e)  and  Eq.  (28),  —  (a,  s) 

for  desired  stations  n  at  each  time  s.  n 

This  change  describes  the  recommended  procedure  for  distributing  the 
total  leading  C„  in  a  reasonable  fashion  along  the  span  of  the  lifting 
surface. 
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